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Abstract. We investigate critical properties of a class of number- conserving cellular 
automata (CA) which can be interpreted as deterministic models of traffic flow with 
anticipatory driving. These rules are among the only known CA rules for which the 
shape of the fundamental diagram has been rigorously derived. In addition, their 
fundamental diagrams contain nonlinear segments, as opposed to majority of number- 
conserving CA which exhibit piccewisc-lincar diagrams. We found that the nature 
of singularities in the fundamental diagram of these rules is the same as for rules 
with pieccwisc-linear diagrams. The current converges toward its equilibrium value 
like £ -1 / 2 , and the critical exponent f3 is equal to f . This supports the conjecture of 
universal behavior at singularities in number-conserving rules. We discuss properties 
of phase transitions occurring at singularities as well as properties of the intermediate 
phase. 
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1. Introduction 

Cellular automata (CA) can be viewed as cells in a regular lattice updated synchronously 
according to a local interaction rule, where the state of each cell is restricted to a finite 
set of allowed values. An interesting subclass of CA consists of rules possessing an 
additive invariant. The simplest of such invariants is the total number of sites in a 
particular state. CA with such invariant, often called "conservative CA" or "number- 
conserving CA" , can be viewed as a system of interacting and moving particles, where 
in the case of a binary rule, l's represent sites occupied by particles, and O's represent 
empty sites. 

In a finite system, the flux or current of particles in equilibrium depends only 
on their density, which is invariant. The graph of the current as a function of density 
characterizes many features of the flow, and is therefore called the fundamental diagram. 
For a majority of conservative CA rules, fundamental diagrams are piecewise-linear, 
usually possessing one or more "sharp corners" or singularities. There exist a strong 
evidence of universal behavior at singularities, as reported in [I]. 

Conservative CA appear in various applications, and some special cases have been 
studied extensively. Rule 184, which is a discrete version of the totally asymmetric 
exclusion process, is an example of such special case El HI El El El E3 EDj- Although 
much is known about this particular rule, dynamics of other conservative rules exhibits 
many features which are not fully understood, and more general results are just starting 
to appear (TTJ E21 ESI EU EH] • In particular, there exists no general result explaining the 
shape of fundamental diagrams for conservative Ca, although promising results have 
been obtained for some special cases by Kujala and Lukka [TH], who studied a class 
of deterministic traffic rules introduced by jT2j. Using their results, we will study 
numerically behavior of singularities for number-conserving rules with fundamental 
diagrams different that previously reported. 

It should be emphasized at this point that all CA rules considered in this paper are 
strictly deterministic. Nevertheless, they are strongly related to stochastic CA models 
of road traffic flow, which have gained widespread attention in recent years [T5| IT§j. 

2. Number-conserving cellular automata 

In what follows, we will assume that the dynamics takes place on one-dimensional lattice 
of length L with periodic boundary conditions. Let Si(t) denote the state of the lattice 
site i at time t, where i £ {0, 1, . . . , L — 1}, t £ N. All operations on spatial indices i 
are assumed to be modulo L. We will further assume that Si(t) £ {0, 1}, and we will 
say that the site % is occupied (empty) at time t if Sj(t) = 1 (sj(t) = 0). 

Let I and r be two integers such that I < < r, and let n — r — I + 1. The set 
{s i+ i(t), Si + i + i(t), . . . , Si +r (t)} will be called the neighbourhood of the site Si(t). Let / 
be a function / : {0, 1}™ — > {0, 1}, also called a local function The update rule for the 
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cellular automaton is given by 

Si(t + 1) = f(s i+ i(t), s i+l+ i(t), . . . , s i+r (t)). (1) 
The CA rule given by will be called number conserving if 



i-l L-l 



f( X i+h x i+l+l, • • • ; ^i+r) — X i (2) 

for all {xqX\, . . . , G {0, 1} L . Note that the above condition simply states that the 
number of occupied sites is constant in time. 

In [211], Hattori and Takesue demonstrated that CA rules are number-conserving if 
and only if a discrete version of a standard current conservation law dp/dt = —dj/dx 
is satisfied. More precisely, CA rule / is number- conserving if and only if for all 
{xi, X2, ■ ■ ■ , x n } e {0, l} n it satisfies 

f(x!,x 2 , ...,X n )- = J(X!,X 2 , ■ ■ ■ - J(x 2 ,x 3 , . . .,x n ), (3) 

where 

n-l -I 

i x 2 1 • • • i x n—l ) -^/(0,0,. ..,0, x ± ,x 2 , ■ ■ ■ ,x n - k ) + y^xj. (4) 

fc=l ' I ' 3=1 

Applying this to all lattice sites we obtain 

f( s i+h • • • 5 s i+r) — Si = J [Si+l, . . . , Sj_)_ r _i) — J • • • , Sj+ r ), (5) 

where we dropped t dependence for clarity. 

The above equation can be interpreted in a similar way as a conservation law in 
a continuous, one dimensional physical system. In such system, let p(x, t) denote the 
density of some material at point x and time t, and let j(x,t) be the current (flux) of 
this material at point x and time t. A conservation law states that the rate of change 
of the total amount of material contained in a fixed domain is equal to the flux of that 
material across the surface of the domain. The differential form of this condition can 
be written as 

?£ = --L (r\ 
dt dx 1 ' 

Interpreting Si(t) as the density, the left hand side of (JHJ) is simply the change of 
density in a single time step, so that (0) is an obvious discrete analog of the current 
conservation law with J playing the role of the current. 

Let us now assume that the initial distribution \i is a Bernoulli distribution, i.e., 
at t — 0, all sites Sj(t) are independently occupied with probability p or empty with 
probability 1 — p, where p G [0,1]. Let us define 

p(i,t)=E ll (s i (t)). (7) 

Since the initial distribution is ^-independent, we expect that p(i, t) also does not depend 
on i, and we will therefore define p(t) = p(i,t). Furthermore, for conservative CA, p(t) 
is t-independent, so we define p = p(t). For the aforementioned Bernoulli distribution 
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we thus obtain p = p. We will refer to p as the density of occupied sites. The expected 
value of the current J(s i+ i(t), s i+ i + i(t), . . . , s i+r _ 1 (t)) will also be ^-independent, so we 
can define the expected current as 

j(p,t) = E^{J{s i+ i{t), s i+/+ i(t), . . . , s i+r _i(t))). (8) 

The graph of the equilibrium current j(p, oo) = lim^oo j(p, t) versus the density 
p is known as the fundamental diagram. It has been numerically demonstrated 
that for conservative deterministic CA the fundamental diagram usually develops a 
singularity as t — > oo, meaning that j(p, oo) is not everywhere differentiable function of 
p. A well-know example is CA rule 184, for which 

s i (t + l) = f(s i - 1 (t),s i (t),s i+1 (t)y (9) 

and / is defined by 

/(0,0,0) = 0, /(0,0,1) = 0, /(0,1,0) = 0, 

/(o,i,i) = i,/(i,o,o) = i, /(i,o,i) = 1, 

/(1,1,0) = 0, /(l, 1, 1) = 1. (10) 
The above definition can also be written in a form (J3J) as 

Si (t + 1) = Si (t) + J(Si-i(t), Si (t)) - J( Si (t), s i+ i(t)), (11) 

where J(xi, X2) = — X2). The graph of the equilibrium current for this rule is shown 
in Figure H^i. The singularity appears at p = 0.5. Quite often, fundamental diagrams 
of number-conserving rules consist of a finite number of linear segments separated by 
similar singularities, as shown in Figure ^D-d. 

Since number-conserving CA rules conserve the number of occupied sites, we can 
label each occupied site (or "particle") with an integer n e Z, such that the closest 
particle to the right of particle n is labeled n + 1. If y n (t) denotes the position of 
particle n at time t, the configuration of the particle system at time t is described by 
the increasing bisequence {2/n(£)}^_oo- W e can then specify how the position of the 
particle at the time step t + 1 depends on positions of the particle and its neighbours 
at the time step t. For example, for rule 184 one obtains 

y n (t + 1) = y n (t) + mm{y n+1 (t) - y n (t) - 1, 1}. (12) 

Equation ([T2|) is sometimes referred to as the motion representation. For arbitrary 
number-conserving CA rule it is possible to obtain the motion representation by 
employing algorithm described in [22] • The motion representation is analogous to 
Lagrange representation of the fluid flow, in which we observe individual particles 
and follow their trajectories. On the other hand, eq. (|TT|) could be called Lagrange 
representation, because it describes the process at a fixed point in space. The Euler- 
Lagrange analogy has been explored in |23j . 

In practice, the choice between Euler and Langrange description of the number- 
conserving CA is usually dictated by convenience. Cellular automata rules which we 
want to consider in this paper are easier to define using the Euler paradigm, as we will 
see in the next section. 
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Figure 1. Examples of fundamental diagrams for number-conserving rules. 
Singularities are denoted by green circles. Diagrams (c) and (d) represent n = 5 
rules with indicated Wolfram's numbers. 



3. Anticipatory driving rules 

In what follows, we will consider a class of CA rules first introduced in ^7j. They 
are defined in terms of two positive integers m (maximum speed) and k (blocking 
parameter). These rules can be interpreted as simplified deterministic models of road 
traffic flow with driver's anticipation. Occupied sites represent cars on a single-lane 
road. All cars move synchronously using the following algorithm. Driver of each car 
first locates the nearest gap (cluster of empty sites) in front of him, the length of this gap 
denoted by g. If the first empty site (i.e., the first site belonging to the aforementioned 
gap) is further than k sites ahead, the car does not move, otherwise, it moves to the 
site i + v, where v = min(g,m). More formally, using Euler's paradigm, and denoting 
position of the n-th car by y n (t), we can write 

y n (t + 1) = y n (t) + [[y b (t) - y n (t) + 1 < k]] mm{g, m}, (13) 

where 

b = min{z G Z : % > n and y i+ i(t) - y^t) > 1}, (14) 
g = y b +i (t) -yb(t) - 1, (15) 
and where [[P]} = 1 if the statement P is true, otherwise [[P]] = 0. 
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Cellular automaton rule defined by (|13fJT5)l will be referred to as as TZ m .k- To 
illustrate dynamics of this rule, let us assume that m — k — 2. Consider, for example, 
the following configuration: 
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The first row represents configuration of cars A, B, C, D at time t, and the second line 
their configuration at time t + 1. Zeros denote empty sites. Car B sees ahead a gap of 
length 2, so it will move by two sites. Car A sees the same gap, so it will move by two 
sites as well. As a result, cars A and B move together as if they formed a single "block". 
This is called "anticipatory driving" - car A anticipates that car B will move. For rule 
lZ m ,k, cars can move in blocks of length up to k. 
Compare this with m = 2, k = 1 case: 
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Car A has not moved, because the first gap ahead of it begins further than k = 1 sites 
away, so A cannot see it, and, as a result, cannot anticipate that B will move. 



4. Known results for lZ m ,i and TZ± : k 

For rules H m ,i, m = 1,2, . . ., the current j(t, p), as defined by eq. ©, can be computed 
explicitly. As shown in [7j, it is given by 

\m + l){t + 1)\ t+i-i/i _ „Wn-i)+i 



kp. t) = i - p - e rhr t + 1 - * r' +1_i(1 _ pr<m,+ *- (16) 

One can then show that in the limit of t — > oo, the equilibrium current j(p, oo) is a 
piecewise linear function of p given by 

j(p ,=o) = (7 ^<v(-+i), (17) 

1 — p otherwise, 

as shown in Figure ^b). 

For p < l/(m + 1) the average velocity of particles at equilibrium is m, i.e., all 
particles are moving to the right with the maximum speed m. The system is said to be 
in the a free-moving phase. When p > l/(m + 1), the speed of some particles is less 
than the maximum speed m. The system is in the so-called jammed phase. 

The transition from the free-moving phase to the jammed phase occurs at p = 
p c = l/(m + 1) called the critical density. At p c , it is possible to obtain an asymptotic 
approximation of (fTBj) by replacing the sum by an integral and using de Moivre-Laplace 
limit theorem, as done in [7j. At p c this procedure yields 



Ul ( (m+l) (m + l)t ' 

j(P^)-j(p,t) = J 2ir{m + i)t (e~™ -e~ - ), 
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and therefore 



j(p, oo) - j(p,t) ~t 



1/2 



(19) 



Here, by f(t) ~ g(t) we mean that lim^oo f(t)/g(t) exists and is different from 0. 

Very similar results can be obtained for 7^1^, k = 1, 2, . . ., if one takes advantage of 
the duality property. Duality in this context means that if cars are moving to the right 
according to the rule TZ m ,k, then empty sites can be treated as another type of particles 
which are moving to the left according to TZk,m- Therefore, if we replace p by 1 — p in 
eq. (jlfij) . we will obtain the equilibrium current for 



As in the case of lZ mj i, rule 72-1^ exhibits two phases represented by linear segments of 
the fundamental diagram, separated by the singularity at p = k/(k + 1). 

5. Rule 1Z 2 ,2 

For arbitrary m and k both greater than 1, explicit expressions for the equilibrium 
current are not known. Nevertheless, Kujala and Lukka ^B] presented an efficient 
algorithm which can determine the steady-state current starting from a given initial 
state. Using that algorithm and the method of generating functions, they obtained 
polynomial equations relating the equilibrium current and the density. 

When m > 1 or k > 1, the fundamental diagram for TZ mt k is different than in the 
case of H m ,i or Hi^- In addition to the free moving phase and the jammed phase, a 
novel phase appears, which we will call an intermediate phase. The smallest values of m 
and k for which one can observe this phenomenon is m = k = 2, and we will use these 
values in subsequent considerations. 

For m = k = 2, the method of fB] method yields 




(20) 



and 




(21) 




2(1 



if p < Pd, 

if p Cl < p< p c2 , 

P) if P > Pc2, 



(22) 



where 
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(23) 
(24) 



Pel 



7 
1 



7 

Pel, 



and where C is a solution of 



16C 2 + p 2 (l - p 2 )(8C 2 - 36C 3 ) + (1 + 27p 2 (l - p 2 ))C 4 - C 5 = 0. (25) 
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Figure El shows the graph of j(p, oo) as a function of p for rule 7^2,2- Two singularities 
at p c i and p C 2 are clearly visible. By singularities we understand values of p for 
which j(p, oo) is non-different iable with respect to p. The two singularities divide the 
fundamental diagram into three parts. When p < p c i, all particles are moving with 
velocity 2, so we will call this a free moving phase. The region p ci < p < p c2 will be 
called an intermediate phase, and the region p > p c2 , in analogy to the fundamental 
diagram described by eq. (fT7|). will be called a jammed phase. 

The existence of singularities in fundamental diagrams of number-conserving rules 
is not new, and their properties have been documented in Q. The most common 
singularity type is a singularity which separates two linear segments of the fundamental 
diagram, just like p = l/(m + 1) for lZ mt i (eq lTTj) . At all such singularities, numerical 
evidence suggests that the relaxation to the equilibrium follows the same power law, i.e. 

j(p,oo)-j(p,t)~i- 1/2 - 

The two singularities in rule TZ 2j2 are different from the singularity observed in 
i, since they separate linear segments of the fundamental diagram from the nonlinear 
segment. Yet surprisingly, they exhibit the same power law behaviour, as we shall see 
in the next section. 

6. Convergence to equilibrium 

In order to define current for 1Z 2j2 , we will first write definition of TZ 2j2 using Lagrange's 
representation. Its local function / appearing in eq. (0) can be written as in a compact 
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Figure 3. Graphs of j(p, oo) — (j(p, t)) as a function of time for rule 7?-2,2- 



form as 

f(Xx, X 2 , X 3l X 4 , X 5 ) =Xi + X 2 X 4 - X X X 3 - X\X 2 X 4 - x 2 x 3 x A + x 2 x 3 x 5 

+xix 2 x 3 x 4 - X2X3X4X5 + x 3 x 4 x 5 , (26) 

with I = —2, r = 2. The local current, as defined by eq. (J3J, becomes 

J(xi,x 2 , x 3 , X4) = Xi + x 2 - x\x 2 X4 - x\x 3 + x\x 2 x 3 x A - x 2 x 3 x A . (27) 

In order to estimate the expected current j(p,t) as defined by eq. (jHJ), we performed a 
series of numerical experiments. We start with a lattice of L = 10 6 sites, where initially 
each site is occupied with probability p and empty with probability 1 — p. After t 
iterations, the average current is then given by 



1 

(j(p,t)) = - ^2j(si(t),s i+1 (t),s i+2 (t),s i+3 (t)). 

i=0 



(28 



Figure El shows graphs of j(p, 00) — (j(p,t)) for three different values of the density 
p: exactly at the singularity (p = p cl ), slightly below (p = p cl — 0.01), and slightly 
above (p = p c \ + 0.01). Due to the symmetry of the fundamental diagram of TZ 2}2 , the 
behaviour in the vicinity of p C2 is the same, therefore we are only considering p ci . We 
can see that at the singularity the power law behavior is observed, as evidenced by the 
straight line in the log-log plot. This is not the case, however, for densities smaller or 
grater than than p c \. For p — p c \, fitting the curve 

j(p,oo)-(j(p,t)) = At a (29) 

to the data set visualized in FigureEl we obtain a = —0.489 ±0.005, which is very close 
to values reported previously for singularities in other number-conserving rules, and very 
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Figure 4. Truncated decay time tt(p) as a function of p for T = 100 (+), T = 200 
(x) and T = 400 (*) for a lattice of 1000 sites. Each point represents average of 100 
runs. 



close to the exact value a = — 1/2 known for TZi tTn . This suggests that singularities in 
7^-2,2 have the same nature as singularities in other number-conserving rules. To confirm 
this observation, we will introduce the notion of the decay time defined as 

oo 

t( p ) = J2j(p^)-j(p^)- (3°) 

t=o 

If the decay is of power-law type, i.e. j(p, oo) — (j(p,t)) ~ t a with a < 1, then the 
infinite sum in eq. (jSOj) should diverge. If, on the other hand, the decay is exponentially 
fast, we should see rapid convergence. For five-input (n = 5) number- conserving CA 
rules with piecewise-linear fundamental diagrams investigated in £Q, the divergence of 
r(p) has been observed only at the singularity, while away from the singularity r(p) it 
quickly converged. Figure H] shows graphs of the truncated decay time defined as 

T 

Mp) = ^2j(p>°°) - j(ft,t) (31) 
t=o 

as a function of p for several values of T. One can see that the decay time diverges at 
critical points p c \ and p C 2, similarly as reported in £Q for other number- conserving rules. 

7. Phase transitions 

Divergence of the decay time at the critical point p c \ is somewhat similar to the critical 
slowing down observed in phase transitions - the closer to the singularity, the longer it 
takes to reach the steady state. In fact, it is possible to define the order parameter for 
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rule 7^-2,2) thus interpreting the singularity at p c \ as the critical point in a second-order 
kinetic phase transition. 

We will first note that the maximum possible current occurs when all particles are 
moving with the maximum speed, which for 7^-2,2 is 2, hence the maximum current is 
equal to 2p. The natural choice of the order parameter would be therefore the difference 
between the maximum possible current and the actual current in the steady statei. . We 
will denote this order parameter by M(p), formally defined as 



The order parameter is zero in the free-moving phase, and becomes non-zero in the 
intermediate phase. 



Using eq. (122(1 . it is possible to obtain series expansion of M(p) around p c \ for 

P > Pel- 



The critical exponent beta, normally defined by M(p) oc (p — p c \Y , is therefore equal to 
1. Again, this is in agreement with the value of j3 obtained for other number-conserving 
rules [T]. 

8. Intermediate phase 

In spite of all similarities to other number- conserving rules, the behavior of 7?.2,2 is 
somewhat unusual, especially in the intermediate phase, where the fundamental diagram 
is nonlinear. 

As mentioned earlier, dynamics of number-conserving CA rules somewhat resembles 
dynamics of the kinematic wave equation, which describes propagation of density waves 
in a continuous medium. For the kinematic wave equation, the slope of the fundamental 
diagram at a given point represents velocity of density waves. It turns out that the slope 
of the fundamental diagram for number-conserving CA can be interpreted in a similar 
way. 

In the free moving phase, the slope of the fundamental diagram given by eq. ()22|) 
is equal to = 2. Compare this to Figure Efa) showing the spatio-temporal diagram 
for p = 0.42, which is in the free-moving phase. If we treat regions of alternating blocks 
11 and 00 as the "background" . . . 0011001100110011 . . ., then white structures can be 
clearly identified in that background. They propagate to the right with velocity 2. These 
are analogs of density waves - in fact, they are regions of local density smaller than 1/2, 
i.e., blocks of zeros longer than 2 or block of ones shorter than 2. Free moving phase 
is dominated by such density waves. On the other hand, in the jammed phase, only 
density waves propagating to the left with velocity —2 remain in the steady state, as 

| Note that this is not an order parameter in the strict sense, i.e., there is no obvious symmetry- breaking 
in the ordered phase. 




(32) 
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Figure 5. Spatio-temporal patterns generated by 7?.2,2 for random initial 
configurations with densities (a) p = 0.42, (b) p = 0.47, (c) p = 0.53, (d) p = 0.58. 
Vertical axis (increasing downward) represents time, while horizontal axis represents 
space. Blacks squares represent occupied sites. 

shown in Figure Efd). This again agrees with the slope of the fundamental diagram 
given by eq. (|22J): in the jammed phase, |^ = —2. 

The intermediate phase is the most interesting one, because its dynamics is not 
dominated by a single type of density waves. Figures Efb) and (c) show two examples of 
spatio-temporal patterns in the intermediate phase. One can clearly see that two types 
of density waves are present. Blocks containing isolated zeros propagate to the right, 
while blocks with isolated ones move to the left. When these two types collide, they 
both are slightly delayed, but after the collision, they continue with the pre-collision 
velocity. One could thus say that these density waves are somewhat similar to solitons 
- they preserve their shape and velocity in a collision. The balance between these two 
types of density waves determines the steady-state current. As it turns out, this balance 
not only depends on the density p, but also on the amount of correlations present in the 
initial configuration. When the initial configuration is described by Bernoulli probability 
measure, i.e., all sites are independently occupied with the same fixed probability, the 
fundamental diagram is given by (j22j) . and the transition to the intermediate phase 
occurs at p cl = | — Yet the location of the singularity p c \ strongly depends on the 
assumption of the initial distribution being Bernoulli. 

To illustrate this, we prepared "clustered initial condition" using the following 
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P 

Figure 6. Order parameter M as a function of p for q = (+), q = 0.2 (x), q = 0.4 
(*), and q = 1 (□). Simulations performed on a lattice of 10 4 sites with Ti-2,2 iterated 
10 4 time steps. Each point represents average of 100 runs. 

algorithm. We start with empty lattice of L sites. In order to produce initial condition 
with pL occupied sites, we repeat the following sequence of steps until the number of 
occupied sites reaches pL: 

(i) select one site randomly among all empty sites and make it occupied by a particle; 

(ii) with probability q, move this particle to the site adjacent to the nearest occupied 
site, and with probability 1 — q, leave it in the initial position. 

In the above, q G [0, 1] is a parameter describing "clustering" of occupied sites. When 
q = 0, we obtain random non-correlated configuration. When q = 1, all occupied sites 
form a single continuous block. 

Figure |S] shows the order parameter M as a function of the density p for four 
different values of q. When q = 0, the critical point is exactly at p c i, as expected. 
However, as q increases, the critical point moves to the right, and when q — 1, it reaches 
0.5. Obviously, the same thing happens to the second singularity, for which one could 
define analogous order parameter. When 5 = 1, these two singularities merge and we 
have just one singularity at p — 0.5. That means that the intermediate phase disappears 
as q — »■ 1. 

This behavior is in sharp contrast with the behavior of number-conserving rules 
with piecewise-linear fundamental diagrams. For these rules, the steady-state current 
depends only on the density, and does not depend on the amount of spatial correlations 
present in the initial condition. 
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We investigated properties of the fundamental diagram of rule 7?.2,2- This rule has 
been used as a representative example of a class of number-conserving rules for which 
the fundamental diagram is known, and for which not all segments of the fundamental 
diagram are linear. We found that the nature of singularities in the fundamental diagram 
of 7^2,2 is the same as for rules with piecewise-linear diagrams. The current converges 
toward its equilibrium value like t -1 / 2 , and the critical exponent (3 is equal to 1. 

These results seem to support a more general conjecture that all singularities in 
number-conserving CA rules exhibit universal behavior. It is interesting to note at 
this point that the critical exponent (3 = 1 is also obtained in equilibrium statistical 
physics in the case of second-order phase transitions with non-negative order parameter 
and above the upper critical dimensionality. This indicates that phase transitions in 
number-conserving CA are all mean-field type, and the observed universal behavior in 
in fact mean-field type behavior. In order to explore this issue further, it would be 
necessary to introduce a field conjugate to the order parameter, and then numerically 
compute other critical exponents such as the exponent 7 characterizing the divergence 
of the susceptibility at the critical point, or 5, characterizing the behavior of the order 
parameter at the critical point when the field approaches 0, similarly as done in [21] for 
the simplified Nagel-Screckenberg model. This problem is currently under investigation, 
and results will be reported elsewhere. 
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